A Novel Path Sampling Method for the Calculation of Rate Constants 
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We derive a novel efficient scheme to measure the rate constant of transitions between stable states 
separated by high free energy barriers in a complex environment within the framework of transition 
path sampling. The method is based on directly and simultaneously measuring the fluxes through 
many phase space interfaces and increases the efficiency with at least a factor of two with respect 
to existing transition path sampling rate constant algorithms. The new algorithm is illustrated on 
the isomerization of a diatomic molecule immersed in a simple fluid. 

PACS numbers: 82.20.Db, 82.20.Sb 
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I. INTRODUCTION 

The calculation of rate constants of activated processes 
dominated by rare events, chemical reactions being a 
prime example, is still one of the major computational 
challenges. As transition rates depend exponentially on 
the activation barrier height, the expectation time for an 
event can exceed current computer capabilities by many 
orders of magnitude. As a result most chemical reactions 
can not be simulated by direct molecular dynamics (MD) 
methods, except those with very low activation energies. 
The conventional way to tackle this time scale problem 
is based on transition state theory (TST) and separates 
the problem in two steps Jl], |^, |J- The first step is 
the calculation of the free energy barrier as function of a 
reaction coordinate, the second stage is the calculation of 
the transmission coefficient by sampling fleeting trajec- 
tories departing from the top of the barrier. If the reac- 
tion coordinate is well chosen, the top of the free energy 
barrier corresponds to points in phase space close to the 
true transition state, and the transmission coefficient will 
have a reasonable value. However, in high dimensional 
complex systems the choice of reaction coordinate can be 
extremely difficult and usually requires detailed a priori 
knowledge of the transition mechanism. Consequently, 
an intuitively chosen but wrong reaction coordinate can 
result in a very low transmission coefficient, and hence a 
statistically inaccurate or immeasurable rate constant. 

Chandler and collaborators (|, ||, [?], ||, ^| devised a 
method for which no prior knowledge of the system is 
needed. This method, called transition path sampling 
(TPS), gathers a collection of trajectories connecting the 
reactant to the product region by employing a Monte 
Carlo (MC) algorithm. The resulting path ensemble 
can be used to elucidate reaction mechanisms, transition 
states and reaction coordinates. The TPS method has 
been successfully used on such diverse systems as cluster 
isomerization, auto-dissociation of water, ion pair disso- 
ciation and on isomerization of a dipeptide, as well a re- 
actions in aqueous solution (see Ref. || for an overview). 
Just as in the conventional case mentioned above, an ad- 
ditional second simulation is needed to determine the rate 
constant within TPS. This simulation combines the path 



sampling method with the umbrella sampling technique 
to estimate the probability to reach the product state 
from the initial reactant state. The final macroscopic 
rate constant is given by aplateau in the time derivative 
of a correlation function . In case of two distinct sta- 
ble states this plateau region should always exist at times 
longer than the typical molecular relaxation time. How- 
ever, when reaction pathways are complex and exhibit 
multiple recrossings, these typical molecular relaxation 
timescales can be relatively long. In that case the TPS 
rate constant calculation is computationally expensive, 
as the path length must exceed these timescales. 

In this paper we improve the efficiency of the TPS rate 
constant calculation on several points by introducing an 
alternative scheme for calculating reaction rates, named 
transition interface sampling (TIS) . The first of these im- 
provements is allowing the path length to vary, so that 
by a well chosen definition of the stable states we can 
limit the length of each path to the strict minimum. Sec- 
ondly, the new method is based on the effective positive 
flux through dividing surfaces or interfaces and is con- 
sequently much less sensitive to multiple recrossings or 
diffusive barrier crossings. Thirdly, the number of differ- 
ent types of Monte Carlo moves is reduced, making the 
implementation of the algorithm conceptually simpler. 

This paper is organized as follows: In Sec. [n] we briefly 
describe the existing algorithms and present the theoret- 
ical derivation for the TIS rate constant expression. The 
implementation of the algorithm is discussed in Sec. III. 
We illustrate the algorithm on a diatomic molecule in a 
fluid of repulsive particles and make a quantitative com- 
parison to the original TPS calculation in Sec. IV. We 
end with concluding remarks in Sec. [y|. 



II. THEORY 

A. Transition state theory and the calculation of 
rate constants 

Consider a dynamical system in which transitions can 
take place between two stable states A and B. If the bar- 
rier between A and B is sufficiently high, the system will 
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show exponential relaxation for which the forward and 
backward rate constants k A B and k^A are well defined 
and can be expressed in terms of microscopic properties. 
Measuring these rate constants by computer simulation 
is traditionally done by the two stage Bennett-Chandler 
(BC) procedure based on the principles of TST J3|, Q. 
The first step is the calculation of the reversible work 
or free energy to bring the system from stable state A 
to the transition state. This free energy ^(A) has to 
be calculated as a function of a suitably chosen reaction 
coordinate A. This A can be a complex function of all 
particle coordinates r and momenta p: A = X(x), with 
x = {r,p}. The maximum in F(X) defines the transition 
state dividing surface A* pTlfl - By convention, the 
system is in A if X(x) < A* and in B if X(x) > A*. 

The main assumption in TST is that any trajectory 
coming from A and crossing the transition state divid- 
ing surface A(x) = A* will remain at the B side of the 
dividing surface for a long time. The reaction rate can 
therefore be expressed as the positive flux through the 
multidimensional dividing surface A* 
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(1) 



where Xt specifies the set of coordinates and momenta 
of the system at time t, the dots denote derivatives with 
respect to time t, the brackets (. . .) denote equilibrium 
ensemble averages and 8(x) and 5(x) are the Heaviside 
step-function and the Dirac delta function respectively. 
In the last equality of Eq. ([[]) the connection to the re- 
versible work F(X) is made, and (3 = 1/ksT, where 
is Boltzmann's constant and T is the temperature. The 
subscript A* to the ensemble brackets, indicates that the 
ensemble is constrained to the top of the barrier A* . 

We consider the system to be completely deterministic 
and thus we can write Xt = f{x t ',t — t') = f(xo,t), in 
which / is the time-propagator function. Evaluation of 
the function f(x,t) requires integrating the equations of 
motion over the time interval t starting with configura- 
tion x. Nevertheless, the equations derived in this paper 
are still valid when applied to stochastic dynamics. 

Even when the TST assumption is accurate, it can be 
extremely difficult to find a proper reaction coordinate 
for which recrossings do not occur. As a result a wrong 
choice for the reaction coordinate will give a much lower 
free energy barrier than the real activation free energy 
and will correspondingly overestimate the rate constant. 
Figure |l| illustrates that Eq. (|l|) overcounts trajectories. 
One can correct for this overcounting by multiplying the 
TST rate constant with the transmission coefficient n(t) 
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FIG. 1: The thin solid curve show the two dimensional free 
energy landscape in contour plot. A is the chosen reaction co- 
ordinate, q represents all other degrees of freedom. A and B 
denote the state regions. The vertical line at A* corresponds 
to the maximum in the free energy function F(X) as is shown 
at the right upper side. The free energy as function of the 
ideal reaction coordinate is also shown at the right lower side. 
This reaction coordinate is a complex function of all degrees 
of freedom Ajdeai = Aideai(A, q) and the corresponding free en- 
ergy function has its maximum at the true transition state 
dividing surface Ajd oa i = A* dcal . This true dividing surface is 
the dashed curved line. The corresponding free energy barrier 
is much more narrow and higher than the artificial barrier due 
to the incorrect reaction coordinate. Four possible trajecto- 
ries are shown. The black solid arrows indicate a positive flux 
through the surface A* and the white solid arrows indicate 
the negative fluxes. TST rate expression (Q) counts all posi- 
tive fluxes of trajectories I, II and III. Consequently, non-true 
reactive events like I and III have a artificial contribution to 
the rate constant and also trajectory II is overcounted one 
time. To correct for this, one can calculate the transmission 
coefficient k. In the TPS equation (^J), if A a = As = A*, 
trajectories III and IV are not counted because of the flA(xo) 
term. Trajectories I and II are correctly counted in the final 
summation due to the cancellation of positive and negative 
flux terms. 



to obtain the true rate constant 

k AB (t) = kl s B T n(t). 



(2) 



The calculation of the time dependent transmission co- 
efficient K,(t) constitutes the second part of the two stage 
BC procedure 0. n{t) belongs to the approximate 
dividing surface A* |l|, Q and can be determined by 
taking an ensemble average of many short trajectories 
starting from the dividing surface: 



«(*) 
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X(x )e(X(x )) 



X(x )9(X(x t ) - X*)) (3) 

/ A* 



After a short molecular time t mo \ the trajectories are 
committed to a stable state and K,(t), and hence kA B {t), 
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become constant: the transmission coefficient k, and the 
rate constant fc^Bj respectively. It is however impor- 
tant to start sufficiently close to the true transition state 
dividing surface. Otherwise the transmission coefficient 
will be extremely low, making an accurate estimate of the 
rate constant problematic or even impossible. In many 
cases, in particular for complex condensed matter sys- 
tems, a sufficiently close reaction coordinate is difficult to 
find and requires considerable a-priori knowledge about 
the system. 



B. Transition path sampling 

Transition path sampling (TPS) is developed to over- 
come the difficulties mentioned above (|[ ^, 0, ||, Its 
main advantage is that no prior knowledge of the tran- 
sition state is needed. The rate constant in TPS is ex- 
pressed as the time derivative of a general time correla- 
tion function. 



We can rewrite the time dependent rate constant of 
Eq. |) into §: 



k T A P B S {t) = ^ t C{t), C(t) 



(h A (xa)h B (x t )) 
(h A {x )) 



(4) 



in which Ha{x) and hs{x) are the characteristic functions 
defined by: 

1ia(x) = 1, if x e A, else Ha{x) = 
fiB(x) = 1, if x e B, else hs(x) = 0. (5) 

In case of a single order parameter X(x t ) describing 
the transition, the phase space regions A and B are 
defined by A^i and Xb- x t € A if \(x t ) < Aa and 
x t £ B if A{x t ) > Xb- Knowledge of the precise loca- 
tion of the dividing surface A*, A^i < A* < Xb, is not 
required in TPS. Therefore, the order parameter A does 
usually not correspond to the reaction coordinate. 

The microscopic expression for the rate constant in 
Eq. (U) is time dependent, while the phenomenological 
rate constant is not. However, just as the transmis- 
sion coefficient n(t) becomes a constant, the time depen- 
dent function k^g S (t) reaches a plateau after a molecu- 
lar timescale t mo \. The phenomenological rate constant 
is equal to the plateau value: k^B — ^^'(T). This 
plateau region should always exist for times T between 
the molecular timescale and the characteristic reaction 
time: t mo \ < T <C t IKn . In other words, T is larger than 
the timescale to commit to one of the stable states, but 
much shorter than the expectation time i lxn of a com- 
pletely new reactive event. If we take A^ = = A* and 
the limit t — > 0+, Eq. (Q) transforms into the expres- 
sion for the positive reactive flux or, equivalently, the 
TST equation ([!]). For t > 0, however, the reactive flux 
measured by Eq. (Q) no longer consists of purely positive 
contributions. The final rate constant is a sum of pos- 
itive and negative fluxes, and thus the overcounting of 
trajectories in Eq. (|l|) is circumvented. (See fig. ([!])). 



L.TPS 
^AB 



(t) 



(h B (t)) 



A,H B (T) 



{hB{t'))A,H B {T) 



C(t% 



(6) 



where H B {T) = max 0<t <T h B {x t ) and (...}a,h b (T) de- 
notes an average on the ensemble of paths of fixed length 
T starting in A and entering B at least once [@. These 
ensemble averages are evaluated using a Monte Carlo pro- 
cedure employing the shooting and shifting moves ||. 
The two factors in Eq. (j|) have to be evaluated sepa- 
rately. First, a path sampling simulation is performed to 
compute {hsit)) a,h b (T) m the interval [0, T]. The path 
length T must be long enough for the time derivative to 
display a plateau. Subsequently, one chooses a t' in inter- 
val [0, T] and computes C(t') using the path sampling in 
combination with an umbrella sampling technique . A 
drawback of the TPS rate constant calculation is that the 
function k^g S (t) can ^ e strongly oscillatory because of 
recrossings and will reach a plateau only after a relatively 
long time. The path length in TPS must exceed the typ- 
ical timescale of these oscillations, and consequently, in 
that case TPS is computationally costly. 



C. Transition Interface Sampling 

Just as the BC and the TPS rate constant algorithms, 
the TIS method is based on a flux calculation. In con- 
trast to these schemes, however, TIS measures the effec- 
tive positive flux |18| , instead of a conditional general flux 
as in Eq. (||) or Eq. ([!]). This implies that only positive 
terms contribute to the rate, allowing for faster numeri- 
cal convergence. A flux is normally defined through a hy- 
persurface in phase space defined by an order parameter, 
the reaction coordinate. But, similar to the TPS case, 
we do not want to suffer from a bad choice of reaction 
coordinate. Therefore, instead of using a single dividing 
surface, we introduce a series of interfaces through which 
we measure this flux. We then derive an expression that 
relates the flux through a certain interface to the flux 
through an interface which is closer to A to replace the 
expensive TPS umbrella sampling procedure. 

In order to formulate a proper flux, we have to divide 
the entire phase space into two complementary regions 
called overall states A and B. These states do not only 
depend on the position at the time of consideration but 
also on its past behavior. Overall state A covers all phase 
space points lying inside stable region A, which consti- 
tutes the largest part, but also all phase space points 
that visit A, before reaching B when the equations of 
motion are integrated backward in time. Similarly, state 
B comprises stable state B and all phase points, coming 
directly from this state in the past, i.e. without having 
been in A. It is useful to generalize the characteristic 
functions in Eq.(^|) for an arbitrary phase space region fi 

hn(x) = 1, iixefl, else h n (x) — 0. (7) 
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For each phase point x and each phase space region f2 we 
can determine the minimum (first entrance) times Vq(x) 
and t^ix) needed to reach f2 starting from configuration 
x by integrating the equations of motion backward and 
forward in time, respectively: 

t b n (x) = -max[{t\hu(f(x,t)) = 1 At < 0}] 
t f Q (x) = + mm[{t\h n (f{x,t)) = lAt>0}], (8) 

where the min and max function return respectively the 
lowest and highest value of their arguments. In addition, 
it is useful to define for each phase point x and each set 
of two non-overlapping phase space regions {f2i, ^2} the 
following characteristic functions: 



1 ifh ni (f(x^t b niun2 (x))) 
otherwise 




_ J 1 if '*ni(/(a ! »+*n 1 un a (*))) = 1 
otherwise 



' (9) 



In words, these functions measure whether a trajectory 
reaches fii before VL 2 or not. As the system is ergodic, 
each phase space region will be visited in finite time and 

thus tyiuth ( x ) + ^Q 2 ,Qx( x ) = h f n u n 2 ( x ) + j&.niO'O = 
1 for any x. Using these definitions the characteristic 

functions for the overall states A and B are given by 



h A (x) = h b AB (x), 



h B (x) = h b B A (x). 



(10) 



These states together span the complete phase space, as 
the system can never stay in the intermediate region be- 
tween A and B forever. The overall states A and B do 
not sensitively depend on the definition of stable state A 
and B as long as it is reasonably. Of course, the stable 
regions should not overlap, each trajectory between the 
stable states must be a true rare event for the reaction 
we are interested in. In addition, the probability that 
after this event the reverse reaction occurs shortly there- 
after must be as unlikely as an entirely new event. In 
other words, the system must be committed to the sta- 
ble states. Therefore, a reasonable definition of A and B 
requires that they should lie completely inside the basin 
of attraction of the respective two states (see also 
Rcf. ||). Special care has to be taken with this condition 
for processes which show many recrossings between state 
A and B before settling down. Such processes can occur 
in solution or in dilute gasses. For instance, for organic 
reactions in aqueous solution, a rare specific hydrogen 
bonded network can lower the bond-breaking barrier and 
initiate the reaction. If the lifetime of those rare solvation 
structures is high, a sudden reverse reaction can occur as 
the barrier for the backward reaction is also lowered by 
the same amount |0, [l4| . A similar phenomenon can 
happen in dilute gasses for which rare spontaneous fluc- 
tuations in the kinetic energy are the main driving force. 
A particle moving from one state to another due to a very 
high kinetic energy as result of sequence of collisions can 
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FIG. 2: Example of phase space regions in TIS. Thin solid 
curves denote the free energy contour lines, qi and qi are two 
arbitrary projections of the degrees of freedom. A and B are 
the two stable states. The dots on the three shown trajectories 
indicate the positions of the system at successive time steps. 
The overall state A and B are indicated by black and white 
dots respectively. Only one trajectory starts in A and ends in 
B and is therefore a true reactive event. The system changes 
from state A into B when it enters region B for the first time. 
It can leave stable region B shortly thereafter, but never go 
back to A in a short time. The stable regions have to be 
chosen to fulfill that condition. 



cross the potential energy barrier several times before it 
will dissipate its energy by a new collision and relax into 
one of the stable states (see e.g. Refs. |l^|). These 
problems can in principle be solved by an adequate choice 
of the stable state definitions. For instance, the definition 
can depend explicitly on the presence of certain hydrogen 
bonds or on kinetic energy terms. 

With our definition of overall states A and B we can 
write down our rate equation in the spirit of Eq. (0) : 



k A i 



h-A(xo)hB(xo) 



(11) 



(h A (xo)) 

where the dot denotes the time derivative taken at t = 0. 
This rate expression does not depend on time although 
the evaluation of the characteristic functions still requires 
integration of the equations of motion. The transition 
from A into B takes place when the system coming from 
A will cross the interface As for the first time (see Fig. |^). 
After this event the system will stay in B. Eq. ( |ll| ) counts 
therefore only the first crossing through interface A B and 
is hence equivalent to the effective positive flux expres- 
sion 



kAB = 



hA(x )X(x )9(X(x ))5(X(x ) - As) 

(h A (x )) 
1 (h A (x )9(X B ~ X(x Q ))0(X(x A t) 



lim , 

At^O At 



(12) 
Ab)) 



(h A (x )} 



Note the similarity with Eq. (|l|). Strictly speaking 
6{Xb — A(iEo)) is redundant in Eq. ( |l2| ) as h A (xo) = 
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if 6{\b — A(a;o)) ^= 1. The last expression in Eq. ( O ) is 
most suitable for a numerical approach with At as the 
time step in a molecular dynamics simulation. Evalua- 
tion of Eq. (|l2| ) requires counting all phase space points 
which at t = are just about to cross interface As in 
one time step and will enter region A before B when 
integrating backward in time starting from xq. Unfortu- 
nately, Eq. ([12]) is not very efficient from a computational 
point of view because only a very small fraction of phase 
points close to interface Xb actually belong to A, leading 
to poor statistics. We can enhance the statistical accu- 
racy by relating the flux through As to the flux through 
an interface closer to A. We therefore introduce a set 
of n non-intersecting interfaces Ai, A2, A3, . . . A„, each in- 
terface Ai closer to A than the next interface A-j+i (see 
Fig. ^) . We define the corresponding phase space regions 
Q\i = { X \M X ) > Ai}. In this way fl\ B is equivalent to 
our stable state B, while £l\ A is the phase space outside 
stable state A. By introducing the following definition 

il™ {xo) ° {Xl ~ X ( x o)) d ( X ( x At) - A,)), (13) 

Eq. ( |l2| ) reduces to 

k AB = ($a,x b ) / (h A ) ■ (14) 

where ($a,A;) denotes the effective positive flux through 
interface Ai . The rate constant is thus equal to the effec- 
tive positive flux through interface A b with the condition 
the trajectories came directly from A. Note again that 
(&a,\a ) / (h-A) is equal to the TST rate expression in 
Eq. (|lj) in case A* = A a = X B - The effective flux ($A,Ai) 
can now be related to the effective flux ("J^a^i) through 
an interface Ai_i closer to A by (see Appendix A) 



(3> A>Xi (x )) 



h L,A( X °U* x($A,A«-i («<>)). 

1 ' *A,A i _ 1 

(15) 

denotes the ensemble average over 



where (. . •)$ A , ; . i _ , 
all phase space points xq for which $A,Ai_i(£o) i= 0. 
The factor ( h' u A{?t ^ ^a,x z _ 1 = is the con- 

ditional probability that a trajectory, coming from A, 
passes Ai, given the fact that it has passed the inter- 
face Ai_i at an earlier time. By recursively substituting 
Eq. (15) into Eq. ( p"4| ) the rate constant can be expressed 
as a product of conditional probabilities: 



k A B 



(hA) 



n 

i=l 
n-1 

]\V{\i +1 \\i)V{\ B \K 



l B,A 



(16) 



P 

B.A 



A.X 



-V(X B \Xi)- 



(h A ) \ B - A /* AlXl ' (h A ) 
This expression is the central equation for TIS. Instead 




FIG. 3: Example of the division of the phase space by inter- 
faces. A and B are the stable state regions with interfaces Xa 
and As. The interfaces Ai . . . A6 correspond to a calculation 
of Eq. ( |l6| ) with n = 6. The dashed lines are the sub inter- 
faces in between. Four trajectories are shown corresponding 
to a "P(A4|A3) ensemble calculation. On each trajectory the 
xo time slice is indicated with a circle. Black circles cor- 
respond to h A q (xo) = and white circles correspond to 
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(x ) 



of just calculating the individual terms in the product 
of Eq. (|l6|) we can equivalently determine a continuous 
crossing probability function "P(A|Ai) for A between Ai 
and Xb ■ This is reminiscent of umbrella sampling where 
a free energy difference is usually estimated as a func- 
tion of a continuous parameter A JhJ . When calculating 
the ensemble average for V(Xi\Xi-i) we can also evalu- 
ate "P(A|Ai_i) for interfaces A between A;_i and A, by 
dividing the phase space into a finer grid of sub inter- 
faces (see fig. (j^)). In this way we acquire useful infor- 
mation without significant extra cost, and, in addition, 
a measure for the convergence of the ensemble averages. 
The final monotonically decreasing crossing probability 
function V(X\Xi) can be obtained by matching the his- 
tograms from the different ensemble simulations. Tech- 
niques commonly applied in umbrella sampling such as 
overlapping windows between two successive ensemble 
averages and the use of biasing functions can also be em- 
ployed here. 



III. THE TRANSITION INTERFACE 
SAMPLING ALGORITHM 

Inspection of Eq. (|l6|) clearly shows that the TIS rate 
constant calculation is also a two step procedure. The 
first step, the effective flux ($a,Ai) / (h A ) can be com- 
puted by simply running a MD simulation starting with 
a configuration in A and counting the number of effec- 
tive crossings. For an interface Ai close enough to stable 
state A one can obtain a statistically accurate value. 

The second part of the calculation consists of evaluat- 
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ing the product of the P(Aj+i|Aj) ensemble averages for 
the different interfaces Aj in Eq. (|T^). Here we need to 
sample all paths from region A to either A or that 
exhibit at least one crossing with interface Xi . The Monte 
Carlo moves in TIS are very similar to the shooting move 
used in the TPS algorithm. The main difference is that 
the backward and forward integration is abandoned as 
soon as the edge of either A or is reached. If the 

new path is accepted there is only one phase point x along 
this path for which ^A.\ii x ) 7^ 0; defining phase space 
point Xq. The shifting moves that were required in the 
original TPS implementation to enable proper sampling 
and improve statistical accuracy are here unnecessary. 

To bootstrap the sampling procedure we first generate 
an initial path that starts in A, then crosses the interface 
Xi and finally ends in either A or (see for more 

details on initial path generation Ref. The phase 

space point xq is then defined as the first crossing point 
of this path with interface A^. Further, let r — int(f/At) 
be the discrete time slice index, and r b = int(t^(xo)/At) 
and = int(t^ un ^ (xo)/At) the forward and back- 
ward terminal time slice indices, respectively. Including 
xq, the initial path then consists of = r b + + 1 
time slices. With these definitions in mind is the TIS 
algorithm as follows: 

1. From the current path with length N^°' choose a 
random time slice r, with — r 6 < r < r* . 

2. Change all momenta of x T At by adding small ran- 
dom displacements 6p from a Gaussian distribu- 
tion. Make sure the total momentum is conserved 



3. In case of a constant energy (NVE) simulation, 
rescale the new momenta to the old energy value 
and continue with step 4. In case of constant tem- 
perature (NVT) accept the new momenta (else re- 
ject the whole TIS move) with a probability flCf] : 



1, exp 



£(41) 



e{4 



(«) x 

At) 



Here, E(x) is the total energy of the system at 
phase space point x. 

4. Integrate equations of motion backward in time by 
reversing the momenta at time slice r, until reach- 
ing either A or Reject in case of else 
continue with the next step. 

5. Integrate from time slice r forward until reaching 
either A or f2> 4+1 . Reject if the entire trial path 
does not cross interface Xi, else continue with the 
next step. 

6. Accept the trial path with a probability 



where is the length of the new path. If ac- 

cepted, replace the old path with the new one. 

7. Reassign xo to be the first crossing point with Xi 
and sample the value of hL a{xq) to measure 
V(X i+1 \Xi). 

8. Repeat from step 1. 

As usual in Monte Carlo schemes, any rejection along 
this route implies counting the old path again in the en- 
semble average. The acceptance probabilities at step || 
and step [| are required to satisfy the detailed balance 
condition (see e.g. Ref ]l0| ). 

Instead of generating a complete path and then accept- 
ing or rejecting accordingly to the probability at step 6, it 
is more efficient to determine a maximum path length in 
advance. Before embarking on the time consuming fourth 
and fifth step, we first take a uniform random number a 
between and I and determine the maximum allowed 
path length by: 



N<& = mt(N^/a). 



(17) 



In this way we can directly stop the integration and re- 
ject the TIS move as soon the path length exceeds 

(n) 

the maximum iVm ax . In the course of the TIS simulation 
the path-length fluctuates. This also means that the av- 
erage path length becomes automatically shorter when 
changing from ensemble average V(X i+ i\Xi) to ensemble 
average V(\i\\i-i) closer to A. 

The algorithm presented here does not require shifting 
moves because there is only one unique Xq phase point 
along each pathway. However, one could consider the 
use of path-reversal moves as they have negligible com- 
putational cost and can sometimes facilitate ergodic sam- 
phng §. 



IV. NUMERICAL RESULTS 

A. The model 

We tested the TIS algorithm on a simple diatomic 
bistable molecule immersed in a fluid of purely repulsive 
particles. Such a system has previously been used in illus- 
trating TPS rate constant calculations JjJ and is therefore 
a good starting point for a comparison between the two 
methods. The system consists of N particles in 2 dimen- 
sions with interactions given by a pairwise Lennard- Jones 
(LJ) potential truncated and shifted at the minimum, of- 
ten referred to as the Weeks- Chandler- Andersen (WCA) 
potential @ 



WCA 



(r) = 



4e[(r/a) 




- 12 - (r/a)- 6 } +e if r < r 
^ if r > r , 

(18) 

where r is the interatomic distance, and r = 2 1 ' 6 a. 
Throughout this section reduced units are used so that 



7 



e and er, respectively the LJ energy and length param- 
eters, as well as the mass of the particles are equal to 
unity. The LJ unit of time (mcr 2 /e) 1 / 2 is therefore also 
unity. In addition, two of the N particles are interacting 
through a double well potential 



Vdw{r) = h 



1 - 



(r-r - w) 



2 1 



(19) 



This function has two minima separated by a barrier of 
height h corresponding to the two stable states of the 
molecule: a compact state for r — ro and extended state 
for r = ro + 2w. For a high enough barrier, transitions 
between the states become rare and the rate constant is 
well defined. Hence, this system provides a useful test 
case for the TPS and TIS algorithms. 

The system is simulated at a constant energy E in 
a simulation square box with periodic boundary condi- 
tions. The total linear momentum is conserved and is set 
zero for all trajectories. The equations of motion are in- 
tegrated using the velocity Verlet algorithm with a time 
step At = 0.002. As in Ref. we focus here on the 
computation of the rate constant for the isomerization 
reaction of the dimer from the compact state to the ex- 
tended state. In the following sectio n we describe general 
simulation details. In Section IV C we discuss the results 



for a system with a high enough barrier to avoid recross- 
ings. Subsequently, we reproduce the simulations from 
Ref. fl7fl in Section IV D. These results do show recross- 



ings, and we discuss the consequences for TPS and TIS. 



computed by means of a straightforward MD simulation 
starting in state A and counting the number of effective 
positive crossings through interface Ai, i.e. when the tra- 
jectory is directly coming from A. The second term in 
Eq. (JLq) is computed using the TIS algorithm of Sec. III. 
The basic requirement is a definition of a set of interfaces 
partitioning the phase space. Between these interfaces 
we defined a finer grid of sub-interfaces to construct the 
crossing probability function 7 5 (A|Ai). As in the TPS 
calculation we adjusted the momentum displacement for 
the shooting move to give an acceptance of about 40%. 

Many parameters are involved in the two methods and 
to compare the relative efficiency we measured the CPU- 
time required for an arbitrary fixed error of 2.5% for each 
step in both the TPS and TIS calculations under the 
same computational conditions (lGhz AMD Athlon). In 
both methods the final rate constant consists of a prod- 
uct of factors which have to be calculated independently. 
For each factor we performed M simulation blocks of N 
Monte Carlo cycles and adjusted N such that after M 
block averages the relative standard deviation of each 
term in Eq. (@) and (JTTJ) was 2.5%. The total CPU-time 
is given by summing the individual 2.5% error CPU-times 
for each factor. The final error in the rate constants is 
obtained by the standard propagation rules using all sim- 
ulation results (i.e. not only the ones for the 2.5% error 
CPU time calculation). 



C. System with High Energy Barrier 



B. Methodology 

The TPS rate constant calculation evaluate s the two 
factors in Eq. (pi) separately as explained in Sec. [IB. The 



first term in Eq. (J^j is the ratio between the plateau value 
of the reactive flux correlation function Qib(T)) a,h b (T) 
and the correction (/i.b(£')) a,h b {t)- The second term 
C{t') requires an umbrella sampling simulation in the 
form of a series of window calculations. An order pa- 
rameter is chosen to define the characteristic functions of 
the stable states and to partition phase space in windows 
for the umbrella sampling. Besides shooting and shifting 
Monte Carlo moves to generate new paths in the tran- 
sition path sampling we also employ a diffusion move 
that shifts the path by one time slice in arbitrary di- 
rection. This move is computationally very cheap but 
increases the statistics of the correlation functions. In 
all our simulations we therefore set the percentages for 
shooting, shifting and diffusion to 5%, 10% and 85%, re- 
spectively. The parameters involved are always gaged 
such that the acceptance ratio is around 40% for shoot- 
ing and shifting moves, ensuring an optimum efficiency 
of the sampling . 

The TIS method involves a direct determination of 
the flux and the calculation of the crossing probability 
functions "P(Ai|Ai_i) between a series of successive inter- 
faces as given by Eq. ([HI). The flux term in Eq. (|l6|) is 



This system had a total number of particles N = 25, 
and a total energy E — 25. The square simulation box 
was adjusted to give a number density of 0.7. The barrier 
height was h — 15 and the width-parameter w — 0.5, so 
that the minima of Vd w (r) were located at r ~ 1.12 and 
r ~ 2.12 while the top of the barrier was at r ~ 1.62. 
In the TPS rate calculation we defined stable states A 
and B as r < ta = 1-5 and r > tb — 1.74, respectively. 
We computed the correlation function (hsit)) a,h b {t) us- 
ing TPS with a fixed path length T = 2.0. ' The cor- 
relation function is shown in Fig. ^ together with its 
time derivative, the reactive flux. The latter function 
clearly displays a plateau. Next, we chose four different 
t' = 0.1,0.3,1.0,2.0 and performed umbrella sampling 
simulations using 8 windows to calculate C(t'). In each 
window we measured the probability to find the path's 
end point r(t') at a certain value of r. These probability 
histograms were rematched and normalized. The final 
probability functions are shown in Fig. ^. Integration of 
the area under the histogram belonging to region B leads 
to C(t') and finally to the rate constant 0. In Table | we 
give the values of the different contributions to the rate 
constant given by Eq. |^, together with the rate constant. 
We report the average relative co mput ation time needed 
to reach the 2.5% error (see Sec. [IV BP in Table |J. 

For the TIS calculations we use the same order param- 
eter r and the same definition for region B, i.e. interface 
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t 



FIG. 4: TPS correlation function (/ib(£)) a,h b (t) (top) and 
its time derivative (bottom) for the system with high energy 
barrier. The error is comparable to line thickness. 




FIG. 5: TPS probability distributions P(r,t') for four t' = 
0.1,0.3, 1.0,2.0 for the high energy barrier. The probability 
P(r, t') is the chance that a path of length t' and starting in A 
will have the end point conformation with a diatomic distance 
r. The graph is the result of the matching of eight window 
calculations. These eight windows are defined as r < 1.19, 
1.18 < r < 1.28, 1.27 < r < 1.35, 1.34 < r < 1.40, 1.39 < r < 
1.47, 1.46 < r < 1.54, 1.53 < r < 1.75, r > 1.74. The errors 
on the histogram points are within the symbol size. 



As is set at r = 1.74. Stable state A was defined by set- 
ting Xa — Ai at r = 1.24. This interface is closer to the 
basin of attraction than the TPS stable state definition 
but yields a higher flux term ($a,Ai } / {hji) an d gives bet- 
ter statistics. Note that the different definition of stable 
state A does not change the final rate constant, as the 
overall state A does not sensitively depends on this def- 
inition. The flux term is calculated by straightforward 
NVE MD. As Xa is equal to Ai every positive crossing 
of this interface is counted in the flux because all trajec- 
tories must by default come directly from A. The con- 
ditional crossing probabilities P(Aj+i|Aj) in Eq. ( |l6|) are 



TPS 



f 
0.1 
0.3 
1.0 
2.0 

TIS 



(h B (T)) AB 

WbWWab 
3300±100 

7.54±0.03 

1.236±0.005 

0.553±0.002 



C*(t')/10~ 13 
0.0018±0.0001 
0.76±0.02 
4.8±0.3 
11.4±0.9 



k A ^ B /w- 13 

6.0±0.5 
5.8±0.1 
5.9±0.4 
6.3±0.5 



($AM)/{hA) P(Ab|Ai)/10" 13 k A -, B /10- 
0.1196±0.0005 49±1 5.9±0.2 



TABLE I: Comparison of rate constants for the high energy 
barrier, computed with TPS at different t' and TIS. Con- 
tributing factors from Eq. (^|) and Eq. (^) are also given. 



calculated for n = 5 interfaces between the stable states 
(see fig. [|). Between these interfaces we impose a finer 
grid to obtain the entire crossing probability function. 
The results for each stage and the final rate constant are 
shown in Table |. The rate constants of both methods 
agree within the statistical accuracy, showing that the 
TIS method is correct. In Table || we give the relative 
computation time to reach the 2.5% error for each term. 
In comparing both methods we have to recall that the 



TPS 



t' 
0.1 
0.3 
1.0 
2.0 
TIS 



(h B (T)) A B 
(h B (t')) AB 

11.0 
0.2 
0.1 
0.1 



Wl W2 W3 W4 W5 W6 W7 W8 Total 

0.01 0.05 0.1 0.04 0.23 0.27 1.3 0.01 13.01 

0.01 0.14 0.28 0.13 0.58 0.43 0.19 0.02 1.98 

1.7 1.7 0.9 0.6 3.0 2.6 6.4 0.2 17.2 

0.03 1.8 4.5 4.4 15.3 8.0 20.3 0.6 55.03 



($A,Ai)/(/iyl) Int Ai Int A 2 Int A 3 Int A 4 Int A 5 Total time 
0.07 0.265 0.09 0.15 0.21 0.215 1 



TABLE II: Comparison of CPU-times required for the 2.5% 
error at each stage for the system with high energy barrier. 
The times are renormalized to the TIS total computation 
time. Wl to W8 denote the different windows used in the 
calculation, Int Ai to Int A5 denote the interface ensemble 
calculations. 
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FIG. 6: TIS crossing probability P(A|Ai) = as 
function of A = r for the system with a high energy bar- 
rier. The function is computed by matching the five inter- 
face ensemble calculations. These interfaces were chosen at: 
Ai = 1.24, A 2 = 1.34, A 3 = 1.40, A 4 = 1.46 and A 5 = 1.52. 
The error on the points is within symbol size. The inset is 
an enlargement in linear scale of the last part of the function. 
We clearly detect a horizontal plateau when approaching As . 



efficiency of TPS depends strongly on the choice of t'. On 
the one hand the umbrella calculation of C(t') is much 
faster for low values of t'. But on the other hand the error 
in the correction term (/is(i')) a.h b (t) increases for lower 
t'. As a result there is an optimum t' for the error/CPU- 
time ratio, in this case approximately at t' = 0.3. Even 
for this optimized situation the TIS calculation is about 
two times faster. One could object that the correlation 
function in Fig. [I] has reached a plateau for t = 1.5 al- 
ready, reducing the TPS computation time by a factor 
3/4. But the choice for a path length T = 1.5 can not 
be taken without a-priori knowledge. The first term in 
Eq. (||) implicitly depends on the path length T. Chang- 
ing T would alter the ensemble and might result in a 
different shape of the flux correlation function. We did 
not check this in detail, but we believe that T cannot be 
chosen much smaller without introducing systematic er- 
rors. Furthermore, we emphasize here that we put much 
more effort in optimizing the TPS algorithm by tuning 
t' , the windows, the ratio between shooting, shifting and 
diffusion moves than we did for TIS. 

Figure shows the histograms of path lengths for each 
TIS ensemble calculation and shows why TIS is faster. 
Sampling paths of fixed length with TPS results in spend- 
ing unnecessary computation time inside the initial and 
final stable regions A and B. In the TIS algorithm instead 
every path is adapted to its minimum length. Bringing 
the interface in closer to A reduces these transition times. 
TIS optimizes itself during the simulation. 



FIG. 7: Histograms P(L) of path length L for each ensemble, 
computed for the system with the high energy barrier. Inset 
(a) is an enlargement of the bottom left area, where windows 
2,3,4 display a second peak. They represent that small frac- 
tion of paths that are able to cross all the interfaces up to the 
rightmost interface and do not have to return to A (cf. the 
trajectories with the white circle in Fig. 3). Inset (b): average 
path length in each window. At variance with TPS the TIS 
algorithm adapts the path length to the ensemble. In going 
from interface 5 to interface 1 one gets closer to state A and 
the path length shortens accordingly. 



D. System with Low Energy Barrier 

In order to compare with previous results, we adopted 
the parameters from Ref . M . The total number of par- 
ticles was N = 9, the total energy was E = 9 and the 
square simulation box was adjusted for a number den- 
sity of 0.6. The barrier height is h = 6 and the width- 
parameter is w = 0.25. Minima are at r ~ 1.12 and 
r ~ 1.62, while the top of the barrier is at r ~ 1.37. 
This barrier is much lower than in the previous section 
resulting in more frequent transitions. An approximate 
rate constant could even be achieved by straightforward 
MD simulations. 

For the TPS calculations we defined the stable states 
A and B by r < r& = 1-30 and r > tb = 1.45, respec- 
tively Qj. Using standard TPS simulation we computed 
the correlation function (hB{t)) a,h b (t) with a total path 
length T = 2 (shown in Fig. [8|). Next, we measured the 
probability histograms to find the paths end point at a 
certain order parameter value r for four different times 
t' = 0.1,0.4,0.8,2.0, using five windows @ (see Fig. |). 
As described in the previous section, matching the prob- 
ability histograms and subsequent integration leads to 
C(t'). The resulting final rate constants, shown in Table 
III, are comparable with the results of Ref. R, but more 
accurate. We will discuss these values after giving the 
results of TIS. 
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FIG. 8: TPS correlation function (/ib(£)).a,h b (t) (top) and 
its time derivative (bottom) for the system with low energy 
barrier. The error is comparable to line thickness. 
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FIG. 9: TPS probability distributions P(r,t') for four t' = 
0.1, 0.3, 1.0, 2.0 for the system with low energy barrier P(r, t') 
as in Fig. ^. The graph is the result of the matching of 
five window calculations. These five window calculations are 
defined as r < 1.22, 1.21 < r < 1.26, 1.25 < r < 1.30, 
1.29 < r < 1.46, r > 1.45 The errors on the histogram points 
are within the symbol size. 



Figure [10| shows that fast recrossings can occur for a 
low barrier, implying that r alone is not sufficient as an 
order parameter to define the stable states in the simu- 
lations. Apparently, this does not effect the TPS results 
much, but it is very important for TIS because of the 
assumption that stable region B is really stable and re- 
crossings do not take place. To ensure the stability of the 
TIS stable states we chose a new order parameter that 
not only depends on the inter-atomic distance r in the 
dimer but also on a kinetic term, given by f. The stable 



TPS 



t' 
0.1 
0.4 
0.8 
2.0 

TIS 



(h B (T)) AB 

47.3±0.2 
2.505±0.007 
1.240±0.003 
0.507±0.001 



C(t')/10" 
1.408±0.007 
2.67±0.01 
5.42±0.05 
13.9±0.2 



k A -, B /w- 

6.67±0.04 
6.68±0.03 
6.72±0.07 
7.03±0.09 



v^K^/io- 5 k A ^ B /w- 

0.2334±0.0003 29.6±0.2 6.90±0.06 



states can then be defined by 

E d (r,r) = ~ + V dw (r) 

x e A if r < 1.37 and E d (r,r) < 1.5 
x e B if r > 1.37 and E d (r, f) < 1.5, 



(20) 



where E d is the sum of the kinetic and potential energy 
of the dimer that has a reduced mass of 1/2. In the {r, f} 
plane these stable states form a D-shape and an inverse 
D-shape regions for A and B respectively (see Fig. fill ). 
Crossing the interface Xa or As implies that the vibra- 
tional energy is decreased below the threshold, E d = 1.5. 
This threshold is made low enough to make fast recross- 
ings to the other state unlikely. However, if we would 



TPS 



f 
0.1 
0.4 
0.8 
2.0 



(h B (T))AE 
{h B (t'))AE 

0.68 
0.4 
0.28 
0.35 
TIS 



Wl W2 W3 W4 W5 Total 

0.03 0.009 0.01 0.1 0.001 0.83 

0.09 0.03 0.04 0.25 0.01 0.82 

0.21 0.07 0.11 1.5 0.04 2.21 

0.28 0.38 0.93 7.27 0.14 9.35 



0.015 



Int Ai Int A 2 Int A 3 Total 
0.085 0.45 0.45 1 



TABLE III: Comparison of rate constants for the low energy 
barrier computed with TPS at different t' and with TIS .in- 
cluding the contributing factors from Eq. (^) and Eq. ([16]), 
respectively. Computation times are reported in units of the 
TIS CPU-time. 



TABLE IV: Comparison of CPU-times required for the 2.5% 
error at each stage for the system with the low energy bar- 
rier. The times are renormalized to the TIS total computation 
time. 
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FIG. 10: Intra-molecular distance of the dimer as function 
of time from a straightforward MD simulation for the system 
with the low energy barrier. Horizontal dashed line at 1.37 
corresponds to the top of the potential barrier. Horizontal 
dashed lines at 1.3 and 1.45 correspond to the TPS state 
definitions of Ref . . Insets are enlargements of four typical 
events on a scale of 10. 1) and 4) correspond to true reactive 
events, A —+ B and B — * A respectively while 2) and 3) are 
non-true, fast recrossing events. In particular, event 3) shows 
capricious behavior with many crossings of the barrier. The 
figure shows a clear separation of timescales, i mo i ~ 1 and 
t rxn ~ 1000. 



have chosen it too low the paths would have become very 
long. We evaluated the crossing probability function in 
Eq. (jlf^) for n = 3 interfaces. The entire crossing prob- 
ability function was obtained by partitioning the phase 
space in sub-interfaces of the form r — A and E e i{r, r) = A 
as shown in Fig. |ll|. Note that in TIS multidimensional 
or multiple order parameters can be used in one simula- 
tion without a problem. This is more difficult in TPS, 
where a proper mapping of the complete phase space is 
required. Figure ^ shows the final rematched crossing 
probability. The monotonically decreasing function tends 
to reach a plateau on approaching the last interface. The 
last two values are not exactly equal but differ by 0.03%, 
indicating that a small fraction of the paths crossing the 
one but last sub-interface still succeed to return to A 
without crossing As. This difference is comparable with 
the chance of a new independent transition (given by the 
rate constant). Note that without the kinetic energy def- 
inition for the stable states Eq. (20), the final crossing 
probability and thus the rate constant would have been 
overestimated by a factor 5/4. 

For the effective flux ($a,Ai) / Cm) calcu lation we per- 

In con- 



formed MD simulations as described in Sec IV B . 
trast to the high barrier case, Ai is not equal A^, and 
not all positive crossings with Ai are effective crossings. 
We counted only the first crossing when the system left 
region A and waited until the system fell back to region 
A before counting a new crossing. As the MD trajectory 
sometimes displayed a spontaneous transition to region 



FIG. 11: One calculated path of the low energy barrier sys- 
tem shown in the {r, f } plane. The vertical solid lines are 
the interface Ai,A2 and A3. The curves A a and As are the 
boundaries of the TIS stable states. The dashed lines are the 
sub-interfaces. The path starts at the dot on Aa and crosses 
the barrier three times before dissipating its energy and re- 
laxing into state B. 

B, we stopped the simulation and started again by re- 
placing the system in a randomized configuration of A. 
Table [II shows the final values and the corresponding er- 
rors of these calculations. The relative computation time 



for each term is detailed in table IV 
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FIG. 12: The crossing probability "P(A|Ai) for the system with 
the low energy barrier. The function is computed by matching 
ensemble calculations with interfaces Ai at r = 1.20, A2 at 
r — 1.26 and A3 at r = 1.32. The inset is an enlargement of 
the final part. The function is converging to a plateau but 
has not yet reached it. The different values of the last points 
are due to the presence of fast recrossings.The error is inside 
the symbol size. 
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FIG. 13: Path length distribution P(L) for each interface 
ensemble in the low energy barrier system The inset shows 
the average path length in each ensemble. 



If we compare the final results of table III we see that 
the efficiency of TIS is more than nine times better than 
the TPS efficiency for t' = 2, and more than two times 
better than TPS value for t' = 0.8. But the TPS t' = 0.1 
and t' = 0.4 efficiencies are about 20 % better than TIS. 
When we compare the rate constants, however, we notice 
that the TPS results for different t' do not agree. Among 
the TPS rate constants only the t' = 2 case is consistent 
with the TIS result. We believe that the t' = 0.1 and 
t' = 0.4 results suffer from systematic errors. For in- 
stance, for the shorter paths the TPS simulations might 
not be completely ergodic. Another explanation might 
be that a path length of T = 2 is too short to allow 
convergence of the reactive flux. In the TIS calculation 
the average path length in the three interface simula- 
tions, from the closest to B to the closest to A, is, re- 
spectively, 7.4, 4.3, and 0.63; much longer than the TPS 
path length (see Fig.[l3]). It is therefore surprising that 
the TPS approach with the simple stable state defini- 
tion and very short paths still gives approximately the 
right rate constant. And indeed, when we computed the 
TPS correlation function with the TIS state definitions 
Eq. (|2"c|), we found that the path length had to be at 
least T = 20 to see a plateau. We think that TPS works 
even with the simple state definitions and the short paths 
because both positive flux and negative flux terms con- 
tribute to Eq. (0). The TPS algorithm collects many 
paths of which some are not real transitions, but fast 
recrossings. The cancellation of positive and negative 
terms of these fast recrossing paths ensure the (almost) 
correct final outcome. In TIS each path must be true 
transition event and contributes as a positive term in rate 
equation (^]), enhancing the convergence. This explains 
that the CPU time for the TIS calculation despite the 



much longer paths is still comparable with TPS one for 
low t'. We note that the path ensemble using the more 
strict stable state definition is of course more useful in 
the analysis of the reaction mechanism. 

For a more accurate comparison of the computation 
time we must keep the systematic errors lower than the 
statistical errors. In other words, we have to make sure 
that the results are converged. To test the convergence 
of the flux correlation function in TPS we can derive the 
following equality from Eq. (^|) : 



(Ww) _ C(t') 

(h B (t"))A,Hs(T) C(t"Y 



(21) 



This equation is valid for any t' ,t" < T if T is large 
enough. We found that the equality does not hold for 
the system with the low barrier, indicating that T is too 
low in the TPS calculation. Further examination of the 
flux correlation function (Hb i^')) a,h b (t) reveals that the 
apparent plateau has in fact a small positive slope. Cal- 
culations for higher values of T suggest that one has to 
increase the path length at least to T = 8 to convergence 
to a plateau. With this in mind we think that the TIS 
computation is about a factor five more efficient than the 
TPS algorithm for the model system with the low barrier. 



V. CONCLUSION 

We developed a novel method, named transition inter- 
face sampling, for the calculation of rate constants based 
on transition path sampling concepts. Just as the origi- 
nal transition path sampling, the new method enables the 
calculation of rate constants of transitions between sta- 
ble states separated by high free energy barriers without 
prior knowledge of the reaction coordinate. The new al- 
gorithm is different in spirit from the rate constant calcu- 
lation that was introduced in Refs||, [jj. In TPS the time 
correlation function (hA(xo)fiB(xt)) / (fiA) is determined 
for a single time using an umbrella sampling scheme fol- 
lowed by a calculation of the reactive flux prefactor in 
a separate path sampling simulation. The path length 
used in this simulation has to be long enough for the 
plateau to be reached. The TIS method advocated here 
calculates the flux correlation directly by measuring the 
fluxes through a number of different interfaces and relat- 
ing the flux through one interface to the next one. The 
big advantage of a flux instead of a correlation function 
is that trajectories going through the interfaces all con- 
tribute to the rate whereas in TPS there are recrossings 
to be counted. In addition, the new method improves 
the original TPS method on a several other points. Once 
the interface is reached, the integration of motion can be 
stopped instead of going all the way to region B. In this 
way the TIS algorithm adapts itself to the optimal path 
length. One does not have to optimize the new method 
as much as TPS, where one has to find the optimal t' 
value and a proper balance between shooting and shift- 
ing. Besides being faster, the concept of calculating a 
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flux comes natural with the rate constant definition, and 
implementation of the algorithm is hence simpler. Also, 
multidimensional or even discrete order parameters can 
easily be implemented in TIS. In the illustrative example 
we showed that we can obtain an increase in efficiency of 
at least a factor of two to five with respect to the TPS 
method used in Ref. M. 

However, one has to be more careful in the definitions 
of the stable states, meaning that stable states have to 
be really stable. In TPS the choice of stable states is 
a bit more flexible as the final rate constant consi sts of 
cancellation of positive and negative terms. In sec IVD 



we showed how this problem for TIS can be solved by 
defining stable regions that explicitly depend on kinetic 
energy terms. 

In summary, we believe that the TIS algorithm can 
make the rate constant calculation of many processes fea- 
sible that were hitherto difficult to obtain. For instance, 
chemical reactions in solution, isomerization of clusters 
and conformational transitions in biomolcculcs. In a fu- 
ture publication we will report on these, more complex, 
applications. 
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APPENDIX A: FLUX RELATION 

In this appendix we show how the effective flux <&A,\i 
can be related to the effective flux $A,Ai_i through an 
interface A^_i closer to A. If at time t = a trajectory 
passes interface Xi while having started in A some time 



earlier, there must always be an unique time when it 
passed interface A^_i for the first time. Therefore we can 
write: 

f tb Aun x . ( x o) 

^A,Xi( x o) = ®A.\i(xo) / dt$ A . Ai _ 1 (a;_ t )(Al) 
Jo 

and hence, 



($a,a 4 (zo)) = 

dt ($A,A 4 _x (*-t)*A,A, (*o)0(*Aun A4 (*<>) - *)) 

dt (^A,x^o)^A^(x t W(t b Aun ^ (x t ) - t)) 



= { $A,A i _ 1 (x ) J di$ A ,A, {x t )6{t AunXi (x t ) - t) 



AuCt x 



(xo) 



dt^A.Xiixt) 



*A,A 4 _ 1 (a:o)^ 4 ,A(a:o))(A2) 



The one but last equation follows because for each phase 
point x and phase space region Q it can be shown that 

t > t f n (x) t b n (f(x,t)) < t e(t b n (f(x,t)) -?j = o 



We rewrite the last expression of Eq. ( A2 ) as a different 
ensemble average: 



^A^-AxoWn A {x ) 



^A,x i - 1 (x )h( lx _ A (x ) 

(xo)) 



X (^A.Ai-iOo)) 



1 L,aM) t * (Qa^-M) , (A3) 

1 ' ®A,\ i _ 1 

where (. . ^ denotes the ensemble average over all 
phase points icq for which ^^Aj-i O^o) =/= 0. The last 



equality gives rise to Eq. ([19 
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